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Abstract 

It is shown that, in the infinite size limit, certain systems of globally coupled phase oscillators 
display low dimensional dynamics. In particular, we derive an explicit finite set of nonlinear ordi- 
nary differential equations for the macroscopic evolution of the systems considered. For example, 
an exact, closed form solution for the nonlinear time evolution of the Kuramoto problem with 
a Lorentzian oscillator frequency distribution function is obtained. Low dimensional behavior is 
also demonstrated for several prototypical extensions of the Kuramoto model, and time-delayed 
coupling is also considered. 
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Because synchronous behavior in large groups consisting of many coupled os- 
cillators has been widely observed in many situations, the behavior of such sys- 
tems has long been of interest. Since the problem is difficult to solve in general, 
much work has been done on the simple paradigmatic case of globally coupled 
phase oscillators. Even in this simple context, however, much remains unclear, 
particularly when considering situations in which a large oscillator population 
interacts with external dynamical systems, or when there are communities of 
interacting oscillators with different community and connection characteristics, 
etc. In this paper we consider an approach that allows the study of the time 
evolving dynamical behavior of these types of systems by an exact reduction to 
a small number of ordinary differential equations. This reduction is achieved by 
considering a restricted class of states. In spite of this restriction, for at least 
one significant example [see preceding article], consideration of our derived or- 
dinary differential equations appears to yield dynamics in precise agreement 
with results obtained from considerations not imposing this restriction. Thus 
we believe that our results may be useful in many other contexts. 



I. INTRODUCTION 



Understanding the generic behavior of systems consisting of large numbers of coupled 
oscillators is of great interest because such systems occur in a wide variety of significant 
applications^]. Examples are the synchronous flashing of groups of fireflies, coordination 
of oscillatory neurons governing circadian rhythms in animals [2 L entrainment in coupled 
oscillatory chemically reacting cellsjsj], Josephson junction circuitsjj], neutrino oscillations 
bubbly fluids [6], etc. A key contribution in this area was the introduction of the following 



model by Kuramoto 



K - 

^(t)M = ^ + -^sin(^(t)-^(t)) , (1) 

3=1 

where the state of oscillator i is given by its phase 8i(t), (i = 1, 2, . . . , N), u>i is the natural 
frequency of oscillator i, and the coupling constant K specifies the strength of the influence 
of one oscillator on another. It has been shown [3, H] that in the N — > oo limit there is a 
continuous phase transition such that, for K below a critical value (K < K c ), no coherent 



behavior of the system occurs (i.e., there is no global correlation between the oscillator 
phases), while above the critical coupling strength (K > K c ), the system displays global 
cooperative behavior (i.e., partial or complete synchronization of the phases). 

Among other problems related to (1) that we shall also consider are the case where there 
is a sinusoidal periodic external drive of strength A added to the righthand side of (1) (see 
Refs.y and [HI), 

JV 

X 



K 

ddi/dt = uoi + — sin (#i - 0») + A sm(nt - 6 t ) , (2) 



and the case where there are several communities of different kinds of oscillators where the 
evolution of the phases 9f(t) of oscillators in community a is given by (see Refs. 



dsr/dt = «r + £ - at £ sin « - n) ■ (3) 

cr' = l j=l 

Here a = 1, 2, . . .,s, X a is the number of oscillators of type a, and K aa i is the strength of 
the coupling from oscillators in community a' to oscillators in community a. For all three 
cases (Eqs. (1), (2), (3)), we are interested in the limit X — > oo. We will also consider such 
problems with time delayed coupling (e.g., 0j(t) — > 8j(t — r) in Eqs. (l)-(3)). 

The problem stated in Eq. (2) was first considered by SakaguchijS]. It can, for example, be 
motivated as a model of circadian rhythm 2| . Circadian rhythm in mammals is governed by 
the suprachiasmatic nucleus that is located in the brain and consists of a large population 
of oscillatory neurons. These neurons presumably couple with each other and are also 
influenced (though the optic nerve) by the daily variation of sunlight (modeled by the term 
in (2) involving A). In lOj, we found numerical and analytical evidence that the bifurcations 
and macroscopic dynamics of (2) with large X appeared to be similar to what might be 
expected for the dynamics of a two dimensional dynamical system. This observation was 
the motivation for the present paper. 



The problem stated in Eq. (3) has been previously considered in Refs. and 12] where 
the linear stability of the incoherent state was investigated along with numerical solutions 
for the nonlinear evolution. 
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II. NATURE OF THE MAIN RESULT 



Considering the limit N — > oo, the state of the oscillator system at time t can be described 
by a continuous distribution function, f(cu, 9, t), in frequency uj and phase 9 for the problems 
in Eqs. (1) and (2) or by f a (uj, 9, t) with a = 1, 2, . . . , s for the problem in Eq. (3), where 



and g{uj) and g a (u) are time independent oscillator frequency distributions. 

Our main result is as follows. For initial distribution functions f(u, 9, 0) satisfying a 
certain set of conditions that we will specify later in this paper, we show that 

(i) the evolution of f(uj,9,t) from f(uj,9,0) continues to satisfy the specified conditions, 

(ii) for appropriate g(oo) [or g a (uj)], the macroscopic system state obeys a finite set of 
nonlinear ordinary differential equations, which we obtain explicitly. 

Concerning (i), we define a distribution function h(ui,9) as a function for which h > 
and d9 J dull = 1, and the distribution functions h(uj,9) satisfying our conditions form 
a manifold M in the space D of all possible distribution functions. What point (i) says is 
that initial states in M C D evolve to other states in M. Thus M is 'invariant' under the 
dynamics. Concerning point (ii), we use the so-called 'order-parameter' description to define 
the macroscopic system state. We define the order parameter (or parameters in the case 
of Eq. (3)) subsequently (Eq. (5)) in terms of an integral over the distribution function / 
(or f a for (3)), where this order-parameter integral globally quantifies the degree to which 
the entire ensemble of oscillators (or ensembles a for (3)) behaves in a coherent manner. 
According to point (ii) the evolution of the order parameters is exactly finite dimensional 
even though the manifold M and the dynamics of the distribution function / as it evolves 
in M are infinite dimensional. 

The macroscopic dynamics we obtain allows for much simplified investigation of the 
systems we study. For example, we obtain an exact closed form solution for the nonlinear 
time evolution of the Kuramoto problem, Eq. (1), for the case of Lorentzian g(u). Our 
formulation will be practically useful if at least some of the macroscopic order-parameter 
attractors and bifurcations of the full dynamics in the space D are replicated in M. In this 
regard, we note that numerical solutions of the system (2) for large iV have been carried 
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out in Ref. lp| . and the resulting macroscopic order-parameter attractors, as well as their 
bifurcations with variation of system parameters, have been fully mapped out. Comparing 
these numerical results for the full system (Eq. (2)) with results for the corresponding low 
dimensional system for the dynamics on M (Eq. (14)), we find that all (not just some) 
of the macroscopic order-parameter attractors and bifurcations of Eq. (2) with Lorentzian 
g{oj) are precisely and quantitatively captured by examination of the dynamics on M. These 
results for the problem given by Eq. (2) suggest that our approach may be useful for other 
situations such as Eqs. (3). Another notable point is that Ref. [lOj also reports numerical 
simulation results for Eq. (2) with large iV for the case of a Gaussian oscillator distribution 
function, g(to) = (27rA 2 ) -1 / 2 exp[— (cu — kj ) 2 /(2A 2 )], and the macroscopic order-parameter 
attractors and bifurcations in this case are found to be the same as those in the Lorentzian 
case (albeit at different parameter values). Thus, at least for problem (2), phenomena for 
Lorentzian g(u>) are not special and should give a useful indication of what can be expected 
for other unimodal distributions giyj). 



III. DERIVATION FOR THE EXAMPLE OF THE KURAMOTO PROBLEM 

We now support points (i) and (ii) for the case of the Kuramoto problem, Eq. (1). 
Following that, we will consider other problems, including those associated with Eqs. (2) 
and (3). Because of its relative simplicity, in this section we use the Kuramoto problem as 
an example, but we emphasize that our interest is primarily in developing a method that 
will be useful in less simple cases, such as the problems stated in Eqs. (2) and (3) (see Sec. 

rin 

IV). Following Kuramotop], |8j, we note that the summation in Eq. (1) can be written as 

where r = N^ 1 ^ exp(i9j). Letting N — > oo in Eq. (1), f(u,9,t) satisfies the following 
initial value problem, 

df/dt + 8/86 {{u + (K/2i)(re~ ld - r*e ie )]f) = , (4) 

/•2tt r+oo 

r= d9 dujfe 19 , (5) 

JO J -oo 

where r(t) is the order parameter, and Eq. (4) is the continuity equation for the conservation 
of the number of oscillators. Note that by its definition (5), r satisfies |r| < 1. Expanding 
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f(u,9,t) in a Fourier series in 9, we have 

{oo 
l + [J2fn(u,t) exp{in9) + c.c 
n=l 

where c.c. stands for complex conjugate. We now consider a restricted class of f n (u,t) such 
that 

f n (u,t) = (a(u,t)) n , 

where |a(u;, £)| < 1 to avoid divergence of the series. Substituting this series expansion 
into Eqs. (4) and (5), we find the remarkable result that this special form of / represents a 
solution to (4) and (5) if 

da/dt + (K/2) (ra 2 - r*) + itoa = , (6) 

+ 0O 

dujat(u), t)g(u) . (7) 

oo 

Thus this special initial condition reduces the ^-dependent system, (4), (5) to a problem 
(6), (7) that is ^-independent. However, we emphasize that Eqs. (6) and (7) still constitute 
an infinite dimensional dynamical system becuase any initial condition is a function of u, 
namely a(u, 0). Performing the summation of the Fourier series using X^nLi xU = ~~ x )-> 
we obtain 

ff n A 9(u) (l-|a|)(l + |a|) 

/(<*>, 9, t) = — — 1 i \ o i a \ i 27T7vi Tvi ' ( 8 ) 



where a = |a|e _4 ^ and r/> real. For \at\ < 1 we can explicitly verify from Eq. (8) that / > 0, 
J d9f = g(u})/2n, and that as \a\ / 1 we have / — > 5(9 — ip)g(uj) / 2ir . In order that Eqs. (6) 
and (7) represent a solution of Eq. (5) for all finite time, we require that, as a(u),t) evolves 
under Eqs. (6) and (7) that |a(u;,t)| < 1 continues to be satisfied. This can be shown by 
substituting a = \a\e~ 1 ^ into Eq. (6), multiplying by e 1 ^, and taking the real part of the 
result, thus obtaining 

d\a\/dt + (K/2)(\a\ 2 - l)Re[re'^] = . (9) 

We see from Eq. (9) that d\a\/dt = at \a\ = 1. Hence a trajectory of (6), starting with an 
initial condition satisfying |a(o;,0)| < 1 cannot cross the unit circle in the complex a-plane, 
and we have |a(u;,t)| < 1 for all finite time, < t < +oo. 

One way to motivate our ansatz, /„ = a n , is to note that the well-known stationary states 
of the Kuramoto model[3, Q], both the incoherent state (/ = g/2ir corresponding to a — 0) 



and the partially synchronized state with |r| = const. > 0, both conform to /„ = a n . Thus 
one view of the ansatz is that it specifies a family of distribution functions that connect 
these two states in a natural way. 

To proceed further, we now introduce another restriction on our assumed form of /: we 
require that a(u,t) be analytically continuable from real u into the complex cu-plane, that 
this continuation has no singularities in the lower half u;-plane, and that |a(u;, t)| — > as 
Im(u) — ► — oo. If these conditions are satisfied for the initial condition, a(u>,0), then they 
are also satisfied for a(u>, t) for oo > t > 0. To see that this is so, we first note that for large 
negative u>i = Im(u), Eq. (6) is approximately da/dt = —\u>i\a, and thus a(u,t) — > as 
LOi — > — oo will continue to be satisfied if a(u, 0) — > as Ui — > — oo. Next we note from[l^] 
that ct(u,t) is analytic in any region of the complex tu-plane for which a(u,0) is analytic 
provided that the solution a(u, t) to Eq. (6) exists. To establish existence for < t < +oo it 
suffices to show that the solution to Eq. (6) cannot become infinite at a finite value of t. This 
can be ruled out by noting that our derivation of (9) with uo now complex carries through 
except that there is now an addition term — Ui\at\ on the left hand side of the equation. Thus 
at \a\ = 1 we have d\a\/dt = u>i\a\ < 0, and we conclude that, if \a(u, 0)| < 1 everywhere in 
the lower half complex tu-plane, then \a(u, t)\ < 1 for all finite time < t < +oo everywhere 
in the lower half complex cu-plane. 

Regarding the initial condition a(u,0), we note that, if |a(a>,0)| < 1 for u real, if the 
continuation a(u, 0) is analytic everywhere in the lower half u;-plane, and if the continuation 
satisfies \a(u, 0) | — > as Ui — > — oo, then the continuation satisfies \a(ou, 0)| < 1 everywhere 



in the lower half complex tu-plane 14j. Examples of possible initial conditions are k exp(-iujc) 
with Re(c) > and \k\ < 1, k/(u — d) with \k\ < Im(d), and / °° k(c) exp(— iuc)dc with 
J~|*(c)|cfc<l. 

We can now specify the invariant manifold M on which our dynamics takes place. It is the 
space of functions of the real variables (u, 9) of the form given by Eq. (8) where \ct(u>, t) \ < 1 
for real a(u,t) can be analytically continued from the real oj-axis into the lower half 
cu-plane; and, when continued into the lower half cu-plane, a(u,t) has no singularities there 
and approaches zero as Ui — > — oo. 

Now taking g(uS) to be Lorentzian 



g(u) = g L (co) = (A/7r)[(w - co o y + A 2 } 



21-1 



7 



we can do the uj integral in Eq. (7) by closing the contour by a large semicircle in the lower 
half d>plane. Writing gi(uj) = (2:/u) _1 [(co> — uj — iA)^ 1 — (uj — uj$ + iA)^ 1 ], we see that 
the integral is given by the residue of the pole at uj = ujq — iA. By a change of variables 
(9,uj) — ► (6 — uj t, (uj — uj )/A), we can, without loss of generality set lu — 0, A = 1. Thus 
we obtain r(t) = a*(—i,t). Putting this result into Eq. (6) and setting uj = —i, we obtain 
the nonlinear evolution of the order parameter r = pe~ % ^ (p > and real): 



dp/dt 



i-\k) p + \ Ki ? 







(10) 



and d(f)/dt = 0. Thus the dynamics is described by the single real nonlinear, first order, 
ordinary differential equation, Eq. (10). The solution of Eq. (10) is, 

-1/2 

(11) 



P(t) 



R 



1 + 



R 



P(0)) 



where R = |1 - {2/K)\ 1 / 2 . We see that for K < K c = 2, the order parameter goes to 
zero exponentially with increasing time, while for K > 2 it exponentially asymptotes to 
the finite value [1 — (2/K)W 2 , in agreement with the known time-asymptotic results for the 
case g = gi (e.g., see Ref. 8|). Plots of the nonlinear evolution of p(t) are shown in Fig. 1. 
Linearization of Eq. (10) yields an exponential damping rate of [1 — (K/2)] for perturbations 




FIG. 1: The order parameter p = \r\ versus time. 



around p = for K < 2, which becomes unstable for K > K c = 2, at which point the stable 
nonlinear equilibrium at p = a/1 — (2/K) comes into existence. For K > K c linearization of 
(10) around the equilibrium p = 



'1 — (2/K) yields a corresponding perturbation damping 
rate [(K/2) — 1]. For g = gi the latter damping rate can also be obtained from the recent 
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stability analyses of solutions of Eqs. (4) and (5) [lCt Il5j. We emphasize that our solution 
for r(t) obeys two uncoupled first order real ordinary differential equations (Eq. (10) and 
dcj)/dt = 0), while the problem for a(u,t) (Eqs. (6) and (7)) is an infinite dimensional 
dynamical system (i.e., to obtain ct(u, t) we need to specify an initial function of lu, ct(u, 0)). 
This is further reflected by the fact that linearization of Eqs. (6) and (7) about their equilibria 
yields a problem with a continuous spectrum of neutral modes 15|, ll6|] . Thus the microscopic 
dynamics in M of the distribution function is infinite dimensional, while the macroscopic 
dynamics of the order parameter is low dimensional. 



IV. GENERALIZATIONS 



a. Other distributions g{oj) 

So far we have restricted our discussion to the case of the Lorentzian Ql^)- We now 
consider 

g(u)=g,(u) = (V2M(u i + l)- 1 , 

which decreases with increasing uo as uj~~ 4 , in contrast to gii^) which decreases as uo~ 2 . The 
distribution g±(u) has four poles at u — (±1 ± i)/V2. Proceeding as before, we apply the 
residue method to the integral (7) to obtain 

r(t) = ±[(l + i)r 1 (t) + (l-i)r 2 (t)) , 

where 

n, 2 = a*((Tl-^)/v / 2,t) 
and r 12 (t) obey the two coupled nonlinear ordinary differential equations, 

dn, 2 /dt + (K/2)[r*rl 2 - r] + [(1 t i)/v^]ri, 2 = . (12) 

Thus we obtain a system of two first order complex nonlinear differential equations. Indeed, 
the above considerations can be applied to any g(u) that is a rational function of uj (i.e., 
g(oS) = Pi(w) I 'P^ioS) where P\{uj) and Pi{oS) are polynomials in uj). The requirement that 
g{ui) be normalizable (f g{uj)duj = 1) and real puts restrictions on the possible Pi t2 (uj); e.g., 
P2{u>) must have even degree, 2m, and all its roots must come in complex conjugate pairs (it 
cannot have a root on the real uj axis). Such a g{oS) has m poles in the lower half tu-plane, 
and application of our method yields m complex, first order ordinary differential equations 
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for m complex order parameters. For instance, for the example, g{oj) = g^iuj), above, there 
are two poles in Im(u) < 0, namely, uj = (±1 — i)/y/2, and these two poles result in the two 
order parameters r\ and r^- 
b. External driving 

We now consider the Kuramoto problem with an external drive, Eq. (2). Again taking 
the N — > oo limit for the number of oscillators, we obtain 



dt 86^ 



+ ^{Kr + A)e~ ie - j.{Kr + A)*e ie 



, (13) 



with r(t) given by Eq. (5). In writing Eq. (13), we have utilized a change of variables 
9 — > 9 + Qt to remove the e lQt time dependence that would otherwise appear multiplying 
the A terms. Again assuming that g(u) is a Lorentzian with unit width A = 1 peaked at 
to — uj , and proceeding as before, we obtain the following equation for r(t), 

dr/dt + ^{(Kr + A)*r 2 - (Kr + A)} + [1 + i(Q - uo )]r = . (14) 

Equilibria are obtained by setting dr/dt = in Eq. (14). Depending on parameters (K, Q, A), 
there are either one or three such equilibria [h|. Also, depending on parameters, there may 
be an attracting limit cycle. Whether the equilibria are attractors for Eq. (14) depends on 
their stability which can be assessed by linearization around the equilibria. The equilibria 
obtained for Eq. (14) and their stability are the same as obtained by the analysis of the full 
system (13) as performed in Ref . [lo| . Furthermore, the bifurcations and stability of the limit 
cycle are the same as numerically found in Ref. [lO] . Thus, for this problem, it appears that 
the important observable macroscopic dynamics is contained entirely within the invariant 
manifold M. 

c. Communities of oscillators 

Turning now to the problem of coupled communities of Kuramoto systems given by 
Eq. (3), we introduce different Lorentzians for each community, 

g°(u)=Tr- 1 [(u-co a ) 2 + Al)- 1 , 

and proceed as before. We obtain a coupled system of equations for the order parameter 
associated with each community a, 

1 s 

dr a /dt + {-iu a + A a )r a + K ™' t r -' r - ~ r -'] = ' ( 15 ) 



2 

a' = l 
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where a = 1, 2, . . . , s. Thus we obtain s complex coupled differential equations where s 
is the number of communities. We conjecture that, for s large enough [e.g., s > 2 or 3] 
and appropriate parameter values, there may be chaotic attracting solutions for Eq. (15). 
It would be particularly interesting to see whether such solutions in M are also attractors 
for the macroscopic order-parameter behavior of the full system (3), e.g., by comparing 
numerical solutions of Eqs. (3) and (15). 
d. Time-delayed coupling 

In applications time delay in the coupling between dynamical units in a network is often 
present. For example, the propagation speed of signals between units is finite (e.g., along 
axons in a neural network), and there may also be an inherent response time of a unit to 
information that it receives. Thus time delay has been extensively studied in the context of 
networks of coupled systems, and in particular for the case of coupled phase oscillators jl?! 



18l . [19| . It has been found for such systems that time delay can substantially modify the 
dynamics, leading to a much richer variety of behaviors. In the context of Eqs. (l)-(3), for 
example, the response of oscillator i at time t to input from oscillator j is now related to 
the state Oj of oscillator j at time (t — Tji) where is the time delay for this interaction. 
Assuming that all the delay times are the same, = r, independent of i and j, the quantities 
Oj(t) appearing in the summations in Eqs. (l)-(3) must now be replaced by 9j(t — r). Again 
such a generalization can be straightforwardly incorporated into our method. For example, 
for the external drive problem (Eq. (2) and Sec. IVb) we have in place of Eq. (2), 

= Ui + f S sin ^ (t - r) - e *w + A sin t fit _ 6 *w ■ (16) 

Going to a rotating frame, O'^t) = 9i(t) — Qt, to' = uj — Q, Eq. (16) becomes 

N 

<It ~' -Y 



^ = 4 + # f>W* " r ) " %® ~ ^ - Asin^'(t) . (17) 



i=l 

The summation in Eq. (17) is 

A' 



K 

iV 



Im | e -KW+nr] ^ e ^(t-r)| = KIm { e -i[«i(*)+n-] r (i _ r )} . (18) 

Thus, to include delay, it suffices to replace the term [Kr(t) + A] in Eqs. (13) and (14) 
by [i^e _iQr r(t — r) + A]. E.g., making this substitution in Eq. (14) and setting A = 0, 
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= ujq yields the following first order delay-differential equation for the order-parameter of 
the standard Kuramoto model with coupling delay, 

^ - y [e- lu)0T r{t - r) - e iuj ° T r*(t - r)(r(t)) 2 ] + r{t) = , (19) 

which returns Eq. (10) for r — > 0. We note that our reduced descriptions with delay (e.g., 
Eq. (19)) are (in contrast to Eqs. (10), (14) and (15)) now infinite dimensional dynamical 
systems. For small |r|, linearizing Eq. (19) about the incoherent state (r = 0), and setting 
r ~ e st yields a dispersion relation for s, 

s + 1 = (K/2) exp[-(s + iu )r) , (20) 

in agreement with Ref. In addition, steady synchronized states can be found (as in 

Ref. jl9|) by setting r(t) = r e tr>t in Eq. (19) and solving the result, 

iV ~ y [e- i{wo+v)T - r 2 e i(uJO+v)T ] +1 = 0, (21) 

for the real constants rj and r$. Furthermore, through linearization of Eq. (19) about r = 
r e trit , our formulation can be used to study the previously unaddressed problem of assessing 
the stability of the steady synchronized states, Eq. (21). 



e. The Millennium bridge problem, Ref. f2u] 



Another example is that of the observed oscillation of London's Millennium Bridge in- 
duced by the pacing phase entrainment of pedestrians walking across the bridge as modeled 
by Eqs. (52) and (53) of the paper by B. Eckhardt et al. [20j. In that case, assuming a 
Lorentzian distribution of natural pacing frequencies for the pedestrians, one can use the 
method given in our paper to obtain an ordinary differential equation for the mechani- 
cal response of the bridge coupled to another ordinary differential equation for the order 
parameter describing the collective state of the pedestrians. 

V. DISCUSSION AND CONCLUSION 

Low dimensional descriptions of the classical Kuramoto problem (Eq. (1)) have been 



previously attempted. An early such attempt was made by Kuramoto and Nishikawa 21] 
who used a heuristic approach resulting in an integral equation for r(t). On the basis of 
their work they predict that for small |r(0)| the order-parameter r(t) initially grows (decays) 
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exponentially in time for K > K c (K < K c ) (later shown rigorously and quantitatively in 



Ref. 



Crawford 



221 ]. using center manifold theory, obtains (Eq. (108) of Ref. 22J) an 



equation of the form dp/dt = a(K — K c )p + bp 3 + 0(p 5 ) for K near K c . Another work of 
interest is that of Watanabe and Strogatz 



23j who consider the case where all oscillators 
have the same frequency for both finite and infinite N. By use of a nonlinear transformation 
of the phase variables 9i(t), these authors show that the dynamics reduces to a solution of 
three coupled first order ordinary differential equations. Thus, while macroscopic behavior 
of order-parameter dynamics has been previously addressed for the standard Kuramoto 
problem, it has, until now, never been demonstrated fully (e.g., without the restriction of 



22| to small amplitude, or the restriction of [23j to identical frequencies). Our paper does 



this and also demonstrates that our technique can be usefully applied to a host of other 
important related problems. 

Our work also suggests other future lines of study. For example, can any rigorous results 
be obtained relevant to whether our macroscopic order-parameter attractors obtained by 
considering / in the manifold M have general validity 24j? Are there interesting qualita- 
tive differences between the behavior for Lorentzian g{oS) as compared to other monotonic 
symmetric oscillator distribution functions g(oj)l What other systems, in addition to those 
discussed in Sec. IV, can our method be applied to? 

We thank B. R. Hunt, M. Girvan, R. Faghih and J. Platig for very useful interactions. 
This work was supported by the ONR (N00014-07- 1-0734) and by the NSF (PHY0456249). 
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Figure Captions 
Figure 1: The order parameter p — \r\ versus time. 
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